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We study the dynamics of the totally asymmetric exclusion process with open boundaries by 
phenomenological theories complemented by extensive Monte-Carlo simulations. Upon combining 
domain wall theory with a kinetic approach known as Boltzmann-Langevin theory we are able to 
give a complete qualitative picture of the dynamics in the low and high density regime and at 
the corresponding phase boundary. At the coexistence line between high and low density phases we 
observe a time scale separation between local density fluctuations and collective domain wall motion, 
which are well accounted for by the Boltzmann-Langevin and domain wall theory, respectively. We 
present Monte-Carlo data for the correlation functions and power spectra in the full parameter range 
of the model. 
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I. INTRODUCTION 

One-dimensional driven lattice gases are an interesting 
field of non-equilibrium statistical mechanics where col- 
lective effects give rise to unexpected non-trivial behav- 
ior such as phase transitions, pattern formation, long- 
range order and anomalous diffusion. In this class of 
problems the totally asymmetric simple exclusion pro- 
cess (TASEP) represents an exactly solvable case which 
proved to model some generic features of several systems 
from a rather diverse range of fields: biological transport 
phenomena (e.g. ribosomes moving on mRNA tracks 
and molecular motors moving along microtubules 0, T 
traffic single file diffusion p| and even economics 

In the TASEP one considers a system of identical 
particles moving uni-directionally with a constant rate 
along a finite one-dimensional lattice with sites labeled 
by i — 1, N (see Fig.^l. The lattice spacing is a = L/N 
with L being the total length of the system. The micro- 
scopic state of the system is characterized by occupation 
numbers rii which are binary variables with only two pos- 
sible values rii £ {0, 1}, i.e. we impose a hard-core re- 
pulsion between the particles. Of particular interest are 
systems with open boundaries, where particles enter the 
system at the left end with a rate a and leave at the right 
end with a rate p. Moreover we use sequential dynamics, 
appropriate for biological systems, where each particle 
moves according to an "internal clock" ; a parallel update 
would be more realistic for vehicular traffic. 

The stationary state (current and density profiles) of 
this driven lattice gas model has been studied in great 
detail over the last years (for recent reviews see e.g. 
Refs. 0) 13)- In open systems it has been found by 
some exact methods 0, 0, that there are three dif- 
ferent phases: a high density, a low density and a maxi- 
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FIG. 1: (a) Cartoon of the TASEP model and (b) phase 
diagram as a function of the entrance and exit rates, a and 
j3, respectively, showing the low density (LD), high density 
(HD) and maximal current (MC) phase. 



mal current phase; compare Fig. ^ The coexistence line 
a = j3 < 1/2 marks a discontinuous transition between 
the high and low density phase, while along the lines 
a = 1/2 < j(3 and = 1/2 < a the transitions are contin- 
uous. There are interesting correlation effects reflected in 
the shape of the density profiles at the boundaries. It is 
only along the disorder line, a + /3 = 1, where the density 
profile is flat. 

The dynamic properties of the TASEP model are 
much less studied. For closed systems exact analyti- 
cal results [H El Q E3 have been derived for the 
largest relaxation time r. It is found that r diverges 
with system size as r oc L z with a dynamic exponent 



2 



of z — 3/2. For open systems no exact results are 
available. The spectrum of relaxation times has been 
studied for small systems using exact enumeration tech- 
niques pH IT?) . More recently density matrix renormal- 
ization group studies have shown that the largest 
relaxation times are finite unless a = f3 = 1/2, where 
r oc L 3 / 2 as for the periodic system; the results also in- 
dicate that z = 3/2 in the whole maximal current phase. 
These studies also confirm the results of a phenomenolog- 
ical approach known as domain wall theory |l7l Il9l | . In 
this coarse-grained description of the dynamics it is as- 
sumed that each particle reservoir at the system's bound- 
ary independently fixes a density. The two domains are 
then joined in the bulk by a phase boundary (domain 
wall) , which performs a random walk due to the random- 
ness of particle flux at the boundaries; the domain wall 
(DW) moves left whenever a particle enters the system 
and right whenever a particle exits the system. Along 
the coexistence line the DW theory accounts even for the 
power spectrum at small frequencies (long time) 20] . 

DW theory does not describe the dynamics of local 
density fluctuations. This is evident for periodic systems 
where the average density profile is flat. In this case it 
has been shown that the dynamics of the TASEP model 
maps onto the noisy Burgers equation |2l| . It is actu- 
ally this mapping which allows for the implementation 
of some exact methods such as the Bethe ansatz [l3j . 
Further prog ress has also been made upon using mode 
coupling |22j ] and renormalization group psl ] theories. 

This correspondence between the lattice gas model and 
a Langevin equation for the periodic case suggests to look 
for a similar mapping for open systems. In this context, it 
is useful to recall a successful method in quantum many- 
body systems, known as Boltzmann-Langevin (BL) equa- 
tions. This approach which was first introduced to de- 
scribe electron transport in the presence of disorder and 
phonon scattering [24,l25l |. describes not only the average 
of the electronic distribution function, but also its fluc- 
tuations. This is achieved by amending the Boltzmann 
kinetic equation with a Langevin source which takes the 
stochastic nature of collisions into account. In the con- 
text of electronic transport through nanostructures, this 
method has recently been widely employed in studies of 
shot noise (for a review see 26] ). A close analogy with 
the lattice-gas model studied in this paper arises for the 
following reason: the requirement of no double occupa- 
tion of sites has its direct analog in the Pauli principle 
which forbids double occupation of electronic states. 

The outline of the paper is the following. In the follow- 
ing section we will discuss how the Boltzmann-Langevin 
approach can be applied for the TASEP model and calcu- 
late the correlation function of the linearized version. In- 
terestingly, it will turn out that the resulting Boltzmann- 
Langevin equation reduces to the noisy Burgers equation 
for a = j3 = 1/2. In section lTTTl we give a short description 
of the Monte-Carlo methods used to analyze the corre- 
lation functions of the TASEP model. Results obtained 
from these simulations are compared with the analytical 



results from the Boltzmann-Langevin and the DW theo- 
ries in section HV1 Finally, we present a short summary 
and some outlook. In Appendix^we give a thorough dis- 
cussion of the behavior of the correlation function right 
at the critical point. 

II. STOCHASTIC EQUATIONS OF MOTION 

At a given time t the microscopic state of the system 
is characterized in terms of occupation numbers {rii — 
0, 1}. The dynamics, that we described as a set of rules 
in the previous section, can be formulated in terms of 
a quantum Hamiltonian representation |27l l28l 129^ . In 
the bulk the corresponding Heisenberg equations for the 
occupation number operators rii(t) have the form of a 
lattice continuity equation 

d t rii(t) = J<_i(f)- Ji(t) (la) 

with the current operator 

Ji(t) = n,-(t) (l-n i+ i(t)) . (lb) 

The effect of the entrance and exit rates is equivalent 
to constant particle reservoirs of density a and 1 — /3 at 
auxiliary sites i = and i = N + 1, respectively. 

There are several levels of approximation in dealing 
with the dynamics of the system. If correlation effects are 
neglected altogether one arrives at a set of rate equations 
for the average particle density Piit) = (n,(t)), which 
have a form identical to Eq. with rii replaced by p^. 
To arrive at these equations one has to take the average 
of Eq. and neglect correlations in the spirit of a mean- 
field or a random phase approximation 

(ni(t)n i+1 (t)) -» (m(t)) (n l+1 {t)) . (2) 

Then, in the stationary limit the rate equations are equiv- 
alent to a nonlinear map 

Pt( 1 ~Pi+i) = J (3) 

with a constant stationary current J. Upon exploiting 
the properties of this map one can easily reproduce the 
full phase diagram of the TASEP 0. Actually, it turns 
out that the phase diagram 9] obtained in this way is 
identical to the one obtained from an exact solution of 
the TASEP in the stationary limit |l(J E3- The den- 
sity profiles obtained from such a mean-field approach 
miss correlation effects, especially in the maximal cur- 
rent phase, and the fluctuations of the domain walls. 

A. The Boltzmann-Langevin approach 

To go beyond rate equations we follow a line of ar- 
guments which leads to what is known as Boltzmann- 
Langevin (BL) equations in studies of non-equilibrium 
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transport in electron systems |24l l25j . The right hand 
side of Eq. (|la|) has a form similar to the collision in- 
tegral for impurity scattering in the Boltzmann equation 
balancing in-going and out-going currents. In order to ac- 
count for fluctuation effects around the stationary state, 
we express both the current and the density as the sum 
of a deterministic and a fluctuating part 



n t w p. + Spi = pi , 
Ji ~ Pi(l - pi+x) + SJi 



(4a) 
(4b) 



Since we will use the BL approach only for those regions 
in the phase diagram where the stationary density is to a 
good approximation spatially constant we may set fa = ~p. 
This applies for both the high and low density phase, but 
not for the phase boundary a — (3 < i, where in addi- 
tion to density fluctuations on small scales we also have 
domain wall motion on large scales. The latter modes 
are obviously not accounted for in the BL formulation. 
One also has to be cautious in the maximal current phase 
where boundary layer profiles decay only algebraically as 
one moves from the boundaries towards the bulk 

Upon inserting Eq. |@J into the equations of motion, 
Eq. QJ, we nn d a coupled set of Langevin equations for 
the density fluctuations at each site of the lattice 



d t 8pi(t) = (1 - p) [Spi-i - Spi] - p [8 pi - Spi+i] 
+5pi(Sp i+ i - 5pi-i) - (SJi - <5Ji_i) 



(5) 



In order to close these equations we still need to specify 
the current fluctuations SJi. This can be done by ex- 
ploiting the fact that the occupation numbers are bi nary 
variables, which immediately implies that J 2 = Ji [38j. 
Hence the variance of the current at a particular site is 
given by Var[Ji] = (Ji) (1 — (Ji)). To be consistent with 
the approximations already made, we set (Ji) « p(l — ~p) 
and finally get 



Var[J]=p(l-p)(l-p(l-p)). 



(6) 



Our final assumption is that correlations in the current 
fluctuations are short ranged in space and time such that 
we can write 



(SJ l (t)SJ j {t')) = Var[J] S i;j S(t - t') . 



(7) 



Note that local current fluctuations are due to the fact 
that each particle advances randomly at a given rate (set 
equal to 1), with an exponential distribution of waiting 
times (in the low density limit). 



B. Gradient expansion 

We will now derive a continuous version of the discrete 
BL equations, Eq. JSJ). To this end we set x = ia and 
introduce fields 4>(x,t) — Spi(t) and rj(x,t) — SJi(t) for 
the density and current fluctuations, respectively. Then 
we get to leading order in a gradient expansion 
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dt<t>(x, *) + («- 24>)dx<t> = M - d x r) 



(8) 



where from now on we measure all length scales in units of 
the lattice spacing a. Eq. has previously been derived 
along similar lines in Ref. p3 • The noise correlations are 
given by 



(r)(x, t)r](x',t')) = AS(x - x')5(t - t') 



(9) 



with an amplitude A = p(l — p)[l — p(l — p)]. We have 
also introduced the collective velocity v = 1 — 2/5, which 
happens to coincide with the expression obtained from 
the exact non-equilibrium fluctuation-dissipation theo- 



rem v — dpJ(p) of an infinite lattice gas 
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Note that 



v changes sign at ~p = i where the stationary current 
becomes maximal. 

The convective nonlinearity <pd x 4> in Eq. |(HJ| can be 
read as a "shift" in the collective velocity due to fluctua- 
tions, which we expect to become important for small v, 
i.e. close to the phase boundaries between the low and 
high density phases and the maximal current phase. For 
densities far away from p = 1/2 we will neglect those 
nonlinearities. Then, as will be discussed in the next 
subsection, one can work out all the correlation functions 
explicitly. These will then be used as a guidance for the 
discussion of the Monte-Carlo results in section llVl 

For p = i Eq. (JSJ is identical to the one-dimensional 
Burgers equation [23J , which can be mapped onto the 
Kardar-Parisi-Zhang (KPZ) 31] equation upon introduc- 
ing a new field h via cf> — d x h. The Burgers equation is 
known to give the following scaling form for the correla- 
tion function C(x,t) = (</>(x, t)<f>(0, 0)) p^ . 



G(x 1 t) = x 2x - 2 F(t/x z ). 



(10) 



Here the roughness exponent \ describes the scaling of 
the width of the interface, and the dynamic exponent z 
characterizes the spread in time of disturbances on the 
surface. In the present case the roughness and the dy- 
namic exponent are known to be |23| 
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and z 



(11a) 



This basically means that the system at the critical point 
relaxes superdiffusively i.e. C(0, t) ~ t~ 1 / z with z < 2. 



C. Correlation functions of the linearized 
Boltzmann-Langevin equation 

The linearized Boltzmann-Langevin equation is most 
conveniently analyzed in Fourier space, where it reads 



1 2 

iuj — ivq + —q 



j,uj) = iqn(q,uj) . 



(12) 



From this one can immediately infer for the correlation 
function C(x — x',t — t 1 ) = (cj)(x, t)4>(x' , t')) in Fourier 
space 



Aq 2 



(uj - vqf + ig4 



(13) 
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and direct space 
C(x,t) = 



A 



y/2n\i 



exp 



(x - vtf 
2\tT 



(14) 



Note that C(x,t) is a Gaussian whose center moves with 
a drift velocity v = l — 2p and which broadens diffusively 
starting from a <5-function at t = 0; height H(t) and 
width W(t) are given by 



A 



(15) 
(16) 



The on-site correlation function decays exponentially for 



Hit) = 

W(t) = 2^\t\. 



C(0,t) 



A 



exp 



\t\ 



(17) 



For v = 0, C(0, t) scales as t -1 / 2 as well as C(0,u>) like 
uj- 1 / 2 . Note also that the static limit of the correlation 
function for the linearized theory is 



lim C(g, oj) 
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4v 2 



(18) 



which is identical to the correlation function for a Lan- 
dau theory in a Gaussian approximation usually found in 
equilibrium thermodynamics . This result suggests 
that linearized BL theory can be viewed as the analog 
of the Gaussian approximation for driven lattice gases. 
The form of Eq. H18|l implies a correlation length of l/2v 
which diverges at the critical point p = 1/2. 



III. MONTE-CARLO SIMULATION METHODS 

Of course, linearized Boltzmann-Langevin theory is 
valid only for very low densities (i.e. low values of a). 
To go beyond this low density limit and test the range 
of validity of the linearized BL approach we have per- 
formed extensive Monte-Carlo (MC) simulations. To this 
end we have chosen the random sequential updating al- 
gorithm by Bortz, Kalos and Lebowitz (BKL- or n-fold 
method) |33L l34J| . Since it keeps a list of all sites which are 
possible candidates for a successful update it is (for the 
present case) faster than conventional methods. More- 
over it constitutes a reliable way to simulate real time 
dynamics and achieve an excellent quality in terms of 
data and computational efficiency in both short and long 
time regimes. 

In a first step one generates a random number X G 
[0, 1), that determines which one of the following moves is 
chosen: a particle entering the system, a particle leaving 
the system and particle at site i jumping to the right. 
Then, for a given move a time interval At is chosen from 
an exponential waiting time distribution, where the decay 
time depends on the size of the list. 



In all of our Monte-Carlo runs we started from a config- 
uration generated according to the steady state distribu- 
tion in order to reduce initial transient effects. After equi- 
libration correlation functions were measured and moving 
time averages over O(10 7 ) time windows were performed. 
Average profiles and correlations do not show any differ- 
ences between moving time and ensemble average giving 
explicit proof for ergodicity of the system. 



IV. DYNAMIC CORRELATION FUNCTIONS 

In this section we analyze the correlation function 
C(x,t) of the density fluctuations <f>(x,t) in space and 
time. From now on we take x — as the central site of 
the system such that the system is confined to the inter- 
val [— L/2,L/2}. This is meant to minimize, at least for 
short times, the influence of the system boundaries. We 
assume that the system is in a stationary state at the 
reference time t = 0. 

As expected from the linearized BL equation, the sim- 
ulations show that the correlation function starts from a 
(5-function peaked at the reference site and then moves 
to the right with velocity v and spreads diffusively as 
time progresses; for negative times it moves to the left 
(see Figs. • Note that our simulations confirm that the 
peak of the correlation function moves with exactly the 
collective velocity v = 1 — 2p. 

In the following we are going to discuss the form and 
the time evolution of the correlation function in the a- 
(3 plane. Our simulations have been performed mainly 
along the anti-diagonal a = 1 — (3. This has the advan- 
tage that the exact steady state profile is perfectly flat, 
such that boundary effects are greatly reduced; moreover 
mean field theory is predicts exactly p = a. 

Before we enter the discussion let us have a closer look 
at the characteristic time scales of the system. Correla- 
tion functions decay on a scale (see Eq. 117(1 



Trelax = 2/w 2 



(19) 



but at the same time the maximum of the correlation 
function moves with a velocity v such that it propagates 
the finite length L of the track in a time 



Tprop 



L/2 



(20) 



This implies that one can observe the relaxation of the 
correlation functions only if T prop /T rc i ax = w-j > 1, i.e. 
for rates a and (3 not too close to the phase boundary to 
the maximal current phase and of course for large enough 
systems. 



A. Low density and high density phase 

Due to particle-hole symmetry we restrict our discus- 
sion to the low density regime. The results for the high 



5 



0.002 



0.0015 



H 0.001 



0.0005 



1 1 1 

(n\ t= " 5 




1 1 1 

t=5 


t=-10 f | 




| ! 1=10 


t=-15 iPVl f 




f \f\ t=15 


- t=-20 ,f\j 1 1 




I V*% t=20 " 









0.015 



0.01 



0.005 



-20 -10 





10 20 








i 1 i 

t=-5 




I 1 I 

t=5 


(6) 






t=-10 / | 




f \ 1=10 


t=-15 1 




| \A 1=15 


t= -20 .-/I ! 




1 f\,,-, t=20 




> 











-10 



10 



20 



FIG. 2: Time and space dependent correlation function for a 
system of 50 sites with rates (a) a = 0.01 and j3 = 0.99 and (b) 
a = 0.1 and (3 = 0.9. Averages are taken over 10 million sam- 
ples. For the diluted system there is good agreement between 
the linearized Boltzmann-Langevin theory (dashed line) and 
the Monte-Carlo simulations (solid line). For denser systems 
the agreement is only qualitative. 




FIG. 3: Peak height of the correlation function for a system 
of N = 50 sites and a series of entrance rates a — 0.005, 0.05, 
0.1, 0.5; we have taken j3 = 1 — a. Averages of the MC data are 
taken over 10 millions samples. Upon approaching the critical 
point (a,/3) = (1/2, 1/2) the peak height shows a power law 
behavior with an effective exponent slowly changing from 1/2 
to 2/3, shown as dashed lines in the graph. 



lower height and are broader than the linearized theory 
predicts. Nevertheless, the peak height measured in the 
simulations still shows the i -1 / 2 scaling of the linear the- 
ory for values of a not to close to a = | ; see Fig. 

As can be inferred from Fig. |31 the effective expo- 
nent describing the peak relaxation slowly crosses over 
from 1/2 to 2/3 upon approaching the critical point 
(a, (3) = (1/2, 1/2) along the anti-diagonal of the phase 
diagram. The exponent 2/3 is identical to the inverse of 
the dynamic exponent 1/z of the non-linear BL equation 
(Burgers equation); see Eqs. J3J and (|11|) . 



B. 



Coexistence line (a = f3 < 1/2) 



density phase are obtained upon simply replacing a by 
1-/3. 

Figures [5] show time series of the correlation function 
for entrance rates a = 0.01 and a — 0.1, respectively. 
In both cases the correlation functions start from a 6- 
function peak at t = which then broadens diffusively 
with the width scaling as t 1 / 2 and the height decreasing 
as i -1 / 2 . Correspondingly the maximum propagates to 
the right or left with velocity v = 1 — 2p for ~p < 1/2 
and p > 1/2, respectively (by particle-hole symmetry). 
As noted above this is an exact result valid for all values 
of the entrance and exit rates. In Fig. we have also 
shown results for negative time to highlight the symmetry 
x — » —x and t — > — t which appears in Eq. 1)14)1. 

For low values of a, which corresponds to the low den- 
sity limit, the results from the linearized BL equations 
explain the Monte-Carlo results quantitatively. The the- 
ory still gives the correct qualitative picture for larger 
values of a but shows significant quantitative deviations. 
The actual shapes of the correlation functions have a 



At the coexistence line the density profile is character- 
ized by a fluctuating domain wall separating a low density 
from a high density phase. The dynamics of the domain 
wall can be described as a symmetric random walk with 
reflecting boundary conditions 0, ^| . In addition to the 
domain wall motion, a collective mode, there are still 
stochastic fluctuations of the density in the low and high 
density wings of the domain wall. Both of these modes 
should be visible in a measurement of the density-density 
correlation function. 

Indeed, the profile of the correlation function shows 
two distinct features; see Fig. 0] At t = there is 
a sharp triangle on top of a much broader triangular 
base. The sharp tip is a result of the local density 
(BL) fluctuations and can be explained as follows. Con- 
sider the correlation function of local density fluctuations 
C«^ L (^) — 011 a lattice. Since there are no 

correlations in the bulk of the system for t — it reduces 
to C^ L (0) = ((nf) - (ni) 2 )5 i: j. Finally, upon using that 
the occupation numbers are binary variables and the av- 
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FIG. 4: Time series of the correlation function C(x,t) ver- 
sus x for a = 0.1 and f3 = 0.1; the times are indicated in 
the graph. The MC data for a system with N = 50 sites 
(solid lines) are compared with a hybrid theory combining lo- 
cal density fluctuations described by the Boltzmann-Langevin 
equation and the collective domain wall motion (dashed lines) . 
Averages in the MC data are taken over 10 million samples. 



erage density in the middle of the system is (no) = 
we obtain C^ L = \5ia- The width of the sharp tip is 
actually a finite size effect resulting from the linear inter- 
polation of the data points. The broader triangular base 
is explained below in the context of DW theory. 

The MC simulations show that the sharp tip quickly 
relaxes and broadens, whereas the shape of the triangular 
base evolves on a much larger time scale. 

One can rationalize this behavior upon combining re- 
sults from the DW and BL theory. We start with a dis- 
cussion of the local density fluctuations. One can derive 
the short time dynamics for the correlation function from 



the lattice version of the BL equation, Eq. (JHJ , imposing 
the initial condition (t = 0) = jS^. As explicitly 
shown in Appendix [E] one finds for the on-site correla- 
tion function C^(t) ~ ~ - + 0(t 2 ), while the near- 
est and next-nearest neighbor correlation functions read 
as C^(t) ~ t and C^(t) ~ t 2 . This explains the fast 
relaxation of the central peak. 

In order to understand the broadening of the triangular 
base we first have to recapitulate some key results of the 
DW theory [HIl^. Since one can model the domain wall 
as a symmetric random walker with reflecting boundary 
conditions at both ends of the system, the conditional 
probability of finding the domain wall at site £t at time 
t given that it was at site £o at time t = reads [23] 



12 2 

^(6|£o) =7 + 1 H e-^'cosA^tcosA^o 

i:even 

^ e-^'sinA^sinA^o, (21) 



2 
I 



z:odd 



where A^ = iir/L and T> = a(l — ot)/(l — 2a) is the 
diffusion constant. Note that this diffusion coefficient is 
smaller than the one of the BL fluctuations (which is 
1), T> /~ a < 1, which explains the time-scale separa- 
tion mentioned above. Averages of an observable O are 
understood as integrals over the random variable £t 

(0(t,t')) = Jdtt [<% t ,P st (&)0(t,t')P({; t \&), (22) 

where -P s t(?) is the stationary probability distribution 
function. In the present case it is simply a constant, 
Pst = l/L. 

If one approximates the density profile of the domain 
wall (DW) by a step function, ip(x,t) = a + (l — 2a)8(x — 
£t), the correlation function can easily be calculated as 



C DW (x,x',t-t') = (Mx,t)Mx',t')) - tydfat)) (Mx',t')) 



2(r-2a) J 

— I? — 



y. 



■ COS X t X COS XiX 1 + Y.i:o 



dd 



■ sin XiX sin XiX 1 



(23) 



Here we are mainly interested in the dynamics at 
time scales t < L 2 /T> (so that the system size is large, 
L 2 > £D, and A^ is infinitesimal), where the domain wall 



has not explored the full system yet. Then the sum in 
Eq. (|23|l can be approximated by an integral, and one 
finds 



C DW (x,t) 



2a) 5 




a; I Erf I 



(24) 
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In the limit t — * this exactly reduces to the profile of the 
broad triangular base in Fig0J If one would be allowed 
to just sum the correlation functions obtained from do- 
main wall and local density fluctuations, this would fully 
explain the initial shape of the correlation function. Of 
course, this is not valid rigorously but seems to be a rea- 
sonable approximation. One may argue that the validity 
of the approximation is due to the time and length scale 
separation between the local density fluctuations and the 
collective domain wall motion. 

In this spirit we assume that the total density fluctua- 
tions <!> can be written as a superposition of local density 
and domain wall fluctuations, t) — 4>{x, t) + ip{x, t), 
and that these fluctuations are uncorrelated, (<j)ip) — 
((j)} {tp}. Then the full correlation function can be written 
as a sum 

C(x,t) = C7 BL (x,t) + C7 DW (x,i) (25) 

with C DW (x,t) given by Eq. (2U and the local density 
correlations C Bh (x,t) are obtained either from the con- 
tinuous or the lattice BL equations depending on the time 
scale. Note that the density fluctuations on both wings 
of the domain wall are the same since the average low 
(/?_ = a) and high (p + = 1 — a) density lead to the 
same noise amplitude A. Hence we may describe these 
local density fluctuations by a BL equation with v = 
and A — a(l — a)(l — a(l — a)). As can be inferred 
from Fig. 0] the corresponding analytical results compare 
reasonably well with MC data. 

A convenient way for visualizing the various dynamic 
regimes resulting from domain wall and local density fluc- 
tuations is the power spectrum 

/M^i(|<I>(0,^)| 2 ) , (26) 

where T is the total time of integration. It is obvious 
from Figs.|S]and[|)]that there are three distinct dynamical 
regimes. 

The DW theory, as described above, fully explains the 
low frequency power law regime I {to) ~ oj~ 3 / 2 . As can 
easily be shown from s peci alizing Eq. (|24|l to x = 0, one 
finds /(w) ~ L _1 oj~ 3 / 2 [2(|. The time window where DW 
theory is valid ranges from the hopping time t\ — 1 /T> to 
the time needed to travel a distance comparable to the 
system size tl ~ L 2 /V (note that L is dimcnsionless). 
For larger times one expects finite size effects. In fre- 
quency space this corresponds to the domain [V/L 2 ,T>}. 

For frequencies larger than lj\ > T> the dynamics is 
dominated by local density fluctuations. Those are well 
described within BL theory. Note that contrary to the 
fluctuations of the domain wall, these local density fluc- 
tuations are independent of the system size; see Figs. |S] 
and IH1 For time scales larger than the microscopic hop- 
ping time of an individual particle (which we have set to 
1), one can use the continuum version of the BL theory. 
Hence for u> < 1 one expects I(lu) ~ lo~ x / 2 which agrees 
very well with our MC data; note that when a = /3«0 
the distinction is clear (Fig. [5J while for a = (3 < 1/2 




FIG. 5: Power spectrum for systems of 200, 100 and 50 
sites and rates a = /3 = 0.1. Averages are taken over 256 
samples. The large frequency behavior is dominated by local 
density fluctuations and well described within a BL theory, 
while the small frequency regime is dominated by domain 
wall fluctuations, as a collective mode. The high resolution 
allows for the identification of a dynamic regime due to the 
discrete nature of density fluctuations at very short time. In- 
set: rescaled power spectrum showing the long time (small 
frequency) regime dominated by the DW dynamics. 




FIG. 6: Power spectrum for systems of 800, 400 and 200 
sites and rates a = j3 = 0.4. Averages are taken over 256 
samples. The large frequency behavior is dominated by local 
density fluctuations and well described within a BL theory, 
while the small frequency regime is dominated by domain 
wall fluctuations, as a collective mode. Inset: rescaled power 
spectrum showing the long time regime dominated by the DW 
dynamics. When the time scale of the DW and BL dynam- 
ics are comparable the separation between the two dynamics, 
although still present, is not as sharp as in Fig. |S] 



not only the time scales but also the amplitudes be- 
come comparable since A = a(l — a)(l — a(l — a)) and 
C DW (0,0) - (1 - 2a) 2 /4 (from Eq.|3J. Therefore (see 
Fig. the distinction between BL and DW regime be- 
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comes clearly visible only at very large time (and large 
systems) . At the critical point the amplitude of the DW 
correlation is identical to zero (see Eq. 123(1 and therefore 
the fluctuations are described by BL in its non- linear ver- 
sion. 

For larger frequencies one has to account for lattice 
effects. If one applies the lattice version of the BL theory 
one finds (see appendix |B"|) 



C k (u) = 



2(1 - cos(!f))A 
,2 + (i_cos(£)) : 



(27) 



In order to obtain the power spectrum Ck(u>) has to be 
summed over all modes numbers k. The dominant con- 
tribution for large frequencies are due to wave vectors 
close to the zone boundary, k = L/2 resulting in a power 



spectrum I(u>) 



which is again well confirmed by 



our MC data (Figs. Eland El . 



V. CONCLUSIONS 

In conclusion we have analyzed the dynamics of the 
TASEP model over the whole parameter range of exit 
and entrance rate with emphasis on the behavior in the 
low and high density regime and the corresponding phase 
boundary. It turns out that most of the dynamics can be 
nicely explained in terms of the combined effect of local 
density fluctuations and collective domain wall motion. 
The dynamics of the domain wall is determined by the 
stochasticity in the entrance and exit of particles at the 
system boundaries. Depending on the parameters this 
yields to a random walk with or without drift towards 
the boundaries. For the description of the local den- 
sity fluctuations we have adopted methods from kinetic 
theories for electronic transport, known as Boltzmann- 
Langevin approach. Both, the Boltzmann-Langevin and 
the domain wall approaches are to a large extent phe- 
nomenological and hence limited in the range of appli- 
cability. Hence we have complemented our studies by 
extensive Monte-Carlo simulations of the TASEP model 
using the BKL algorithm which allows us to study the 
real time dynamics with good accuracy. Our main find- 
ings are as follows. For very low densities, the linearized 
Boltzmann-Langevin theory accounts quantitatively for 
the shape of the density-density correlations. It becomes 
less accurate for densities approaching the maximal den- 
sity of 1/2 as expected from the approximate nature of 
the theory. Analogous arguments apply for very high 
densities by virtue of particle-hole symmetry. For densi- 
ties close to 1/2 linearized Boltzmann-Langevin theory is 
quantitatively wrong but still captures the main features 
qualitatively. Exactly at the critical point, a = /3 = 1/2, 
the full Boltzmann-Langevin theory is identical to the 
noisy Burgers equation which is known to be in the same 
universality class as the TASEP model right at this point. 

As summarized by the power spectra in Figs.El& n dEl 
there is a time scale separation between the domain wall 



motion and the local density fluctuations. For frequen- 
cies larger than the hopping time of the domain wall T> 
it is the local density fluctuations which dominate the 
spectrum. Upon using the continuous and the discrete 
version of the linearized Boltzmann-Langevin approach 
we can fully account for the crossover from w -1 / 2 to uj~ 2 
in the spectrum. For low frequencies uj <T> domain wall 
theory gives a power spectrum of lu~ 3 / 2 in agreement 
with the Monte-Carlo data. 

In summary, two rather elementary approaches, do- 
main wall and Boltzmann-Langevin theory, seem to cap- 
ture most of the observed dynamics of the TASEP model. 
This suggests that it may be worthwhile to look for more 
complex systems which also could be described by these 
simple methods. 
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APPENDIX A: CRITICAL POINT a = = 1/2 

For completeness we shortly discuss our results for the 
correlation function right at the critical point (a — (3 — 
1/2). As can be inferred from Fig. 0the temporal evo- 
lution of its shape is qualitatively and quantitatively dif- 
ferent from the low and high density phases. 

The critical exponents are obtained from a finite size 
scaling analysis of the height and width of the correlation 
function. The insets of Figs.Ela and b show Monte-Carlo 
data for system sizes N = 10, N — 25 and N = 50. These 
data can be replotted upon using the finite size scaling 
relations for the height and width, respectively, 

C(x = 0, t) = L 2 ^ 2 g{t/L z ) = L^git/Li) , (Al) 



W(t) 



L**Pf(t/L') = Lf(t/tf), 



(A2) 



which will give us numerical values for the critical expo- 
nents. 

The critical exponents were determined using an al- 
gorithm provided by the authors of Ref.J^|. This code 
computes and minimizes a sum which weights the dis- 
tance from an interpolating function based on all the 
given sequences of data. Errors are extracted measuring 
the width of the minimum of such function (which has 
been tested to be zero if the values are exact). From this 
we obtain for the autocorrelation (peak) 2\ — 0.98 ±0.03 
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FIG. 7: Time series of the correlation function C(x,t) versus 
x at the critical point (a = 0.5 and /3 = 0.5); the times are 
indicated in the graph; the system size is N = 50. Averages 
are taken over 10 million samples. 
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FIG. 8: (a) Width of the correlation function and (b) auto- 
correlation. Averages are computed using 10 7 samples at the 
critical point a = j3 = 1/2. The data are rescaled according 
to Eqs. I|A20 and l|Al|l with the exponents presented in the 
text. Insets show the plots before rescaling. 



and z = 1.52 ± 0.02, and consistently 2\ = 0.92 ± 0.09 
and z = 1.53 ± 0.08 for the width. It constitutes direct 
numerical evidence for the system belonging to the KPZ 
universality class (2% = 1 and z — 3/2), as expected 
from earlier analytical results for periodic systems. Our 
measurements confirm the numerical results in Ref. 1371. 



APPENDIX B: ANALYSIS OF THE 
BOLTZMANN-LANGEVIN EQUATION ON A 
DISCRETE LATTICE 

In this appendix we study the time behavior of the 
correlation function at short times, in the regime where 
the discreteness of the system plays a major role. We 
use 4> n (t) — Sp n (t) as the discrete equivalent of the field 
4>{x, t). We can write an equation of motion for the cor- 
relation function multiplying Eq. JSJ by (/>° : 

= (^0) (1 _ a) _ M) + 

+a (<t> n+ l<t>° n ) - (tint*) + (Vn-l€) (Bl) 

where ~p = a and we neglect the non linear terms. This 
system of equations involves two point correlation func- 
tions for three different lattice sites, but can be written 
in a closed form (at least for short time regimes) by as- 
suming that: 

1. (vnfin) = which makes sense being the noise in- 
dependent on the dynamics itself 

2- (<An-20n) = which is reasonable for short time 

3. (^m^n) ~ \&mn (this originates the sharp tip in 
C(x,t)) 

The central site of the system will be considered the ref- 
erence site (n = 0). Defining C*(i) = (0 n+i (<)</>° ) we 
rewrite Eq. IjBlfl as: 



dC 
dt 



(1 - a)C_i - Co + aC+i 



(B2) 



multiplying Eq. (jBll for 4> n -i and 4>n+i by analogous rea- 
soning we find a system of linear differential equations: 



dv(t) 
dt 



Mv(t) 



(B3) 



with v = (C_i, C 0j C+i)* and M 

For the initial condition vq — 
v(t) = cxp(Mt)v leads to 

1 



-1 a 

1 — a —1 a 

1 - a -1 

(0,1/4,0)* the solution 



C (f) = ^-e^cosh (ty/2a(l - af) 

\-\t+\{\+a-o?)e + 0{t3) (B4) 
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At short time the autocorrelation decays linearly in time 
from a constant value (1/4), while the other terms grow 
linearly: 



c - l{t) = Vur^) 6 '* sinh (y 2 ^ 1 - a 



(B5) 



C+i(i) - C7+i 



1 /l -a 



2a 



sinh ( £y2a(l — a] 



(1 — a) 1 — a 9 „ . q, 
^— -i* - ——t 2 + 0(t 3 ) . 



(B6) 



Note that even relaxing hypothesis (2), assuming there- 
fore C±2 and dealing with a larger matrix M, one 
does not find correction to the leading behavior in time 
for Co(t), since correlation functions for more distant 
sites, as C±%, scale as C± 2 ~ t 2 - 

In order to look at the behavior in frequency space, 
we apply the BL scheme in the discrete lattice and ex- 
trapolate the regime of the correlation function at large 
lj. Let us start from the real space-time Boltzmann 
Langevin Eq. JSJ and let us introduce the discrete Fourier 



transform <fik = Yln=-L/2 <\>n& tknlx ^ L where k indicate the 
mode number. In order to express the linearized BL 
equation in discrete Fourier space, we multiply Eq. (0 
by e jfc ™ 7I 7 i an( j sum over n Even though the system is 
not translational invariant, in this limit the system can 



be considered as infinite and we do not take care of the 
boundaries. Performing a Fourier transform in time we 
get the discrete equivalent of Eq. © : 



i(u> + 2a sin — ) + (1 — e" 

E 



<f>k(w) = %M(1 - e''*) 
(B7) 



and find the correlation function 



C k (cj) 



2A(1 - 



cos ^S) 



l / 

(w-«sin^) 2 + (l-cos^) 2 



(B8) 



where we use the notation v = 1 — 2a and A — a(l — 
a(l — a)) as done above. 

The dominant contribution for large frequencies are 
due to wave vectors close to the zone boundary, k = L/2 
resulting in a power spectrum I (to) ~ oj~ 2 , which is again 
nicely confirmed by our MC data (Figs. [S] and 

The autocorrelation is the sum over all the modes, but 
the dominant contribution for large frequencies are due 
to wave vectors close to the zone boundary, k = L/2: 

1 ^2A(1 -cos fe) , , s 

C(x = 0, w) ~ - L — 2 — = 2Alj ( B9 ) 



fc=0 



which (by Wiener-Khinchin theorem) is the power spec- 
trum mentioned in Sec lIVBl 
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